Decoherence in Weakly Coupled Excitonic Complexes 
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Equations of motion for weakly coupled excitonic complexes are derived. The description allows 
to treat the system in the basis of electronic states localized on individual chromophores, while at 
the same time accounting for experimentally observable delocalization effects in optical spectra. The 
equations are show to be related to the well-known Forster type energy transfer rate equations, but 
unlike Forster equations, they provide a description of the decoherence processes leading to suppres- 
sion of the resonance coupling by bath fluctuations. Linear absorption and two-dimensional photon 
echo correlation spectra are calculated for simple model systems in homogeneous limit demonstrat- 
ing distinct delocalization effect and reduction of the resonance coupling due to interaction with the 
bath. 
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I. INTRODUCTION 

Excitonic interaction determines spectroscopic and 
functional properties of many naturally occurring, as well 
as artificially synthesized macromolecular systems 
The great variability of photosynthetic antennae of plants 
and bacteria that involve only a limited number of dif- 
ferent types of small molecules as building blocks, is to 
a high degree enabled by the large influence of inter- 
chromophore interactions on spectral and energy transfer 
properties of closely packed protein-chromophore com- 
plexes (1, Q. In recent years, some of the properties 
and function of photosynthetic antennae were linked both 
theoretically 0, Q and experimentally [9l-[ll| to various 
types of coherence effects. In some photosynthetic com- 
plexes, dynamics related to long-lasting electronic coher- 
ences was observed during relaxation of excitation energy 
at 77 K and even at room temperature [T3 |. and it 
was speculated that the coherent mode of energy transfer 
can improve the robustness of the energy transfer process 

[13 ■ While quantitative aspect of this improvement is 
a matter of ongoing research, another closely related type 
of coherence, the one accompanying delocalization of ex- 
cited states in weakly coupled chromophore complexes, 
has been found to play a significant role in spectra and 
energy relaxation rates of some bacterial antenae. For the 
peripheral bacterial light harvesting complexes LH2 and 
LH3, it has been shown that event at very weak resonance 
coupling the absorption spectra, relaxation rates and co- 
herent two-dimensional electronic spectra show charac- 
teristic features of delocalization 0, llj- 

Proteins fix the positions of chromophores in a protein- 
chromophore complex, and thus determine their mutual 
interactions. However, the protein surrounding also tune 
local excitation energies of the chromophores to achieve 
further flexibility of the antennae (T5| . The protein en- 
vironment influences the electronic degrees of freedom 
(DOF) also on the ultrafast time scale. The shape of the 
absorption line of a separate molecule is defined by the in- 
teraction of electronic transition with the intramolecular 



vibrations and vibrations (phonons) of a molecular sur- 
rounding. The concept of line-shape function originates 
from pioneering theoretical work by Lax on absorption 
spectrum of a two-level system coupled to harmonic os- 
cillators [l^ . Further development of the theory has been 
based on the stochastic approach [13, EEli description of 
anharmonicity [l9| and the Brownian oscillator model 
for the phonon system [2^. In molecular aggregates the 
intermolecular interactions compete with the electronic 
transition coupling to the intramolecular vibrations and 
phonons, which both are usually considered as thermal 
bath fluctuations. 

Two-limiting cases of relative strength of the resonance 
interaction and the electron-phonon coupling are widely 
used in various theoretical methods. In the strong ex- 
citon coupling regime (with respect to electron phonon 
interaction) the homogeneous bandwidths of the absorp- 
tion/emission spectra and spectral dynamics are related 
to the exciton relaxation (dephasing) caused by the ther- 
mal bath fluctuations. Excitonic splitting of the transi- 
tion energies due to the resonance interaction schemat- 
ically illustrated in Fig. [T] dominates absorption spec- 
trum. In the opposite regime, the intermolecular res- 
onance interaction might be considered perturbatively. 
The excitation dynamics is then characterized by hop- 
ping between the chromophore molecules, and can be well 
described by the Forster theory. Forster resonance en- 
ergy transfer (FRET) is nowadays widely accepted as a 
molecular ruler in various biological systems pll. [2^. 

The above mentioned photosynthetic complex LH2 is a 
good example of a system where both these limiting cases 
persist. This complex is arranged as a highly symmet- 
ric ring of 9 (or 8 depending on the species of bacteria) 
protein-chromophore subunits, each containing two heli- 
cal trans-membrane polypeptides, the a-polypeptide on 
the inner side and the /3-polypeptide on the outer side of 
the ring fsl . The carboxy-terminal domain of this protein 
binds, in the hydrophobic membrane phase, a ring of 18 
(or 16 depending on the species of bacteria) tightly cou- 
pled bacteriochlorophyll (Bchl) molecules with a center- 
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to-center distance of less than 1 nm between neighboring 
chromophores. This ring is responsible for the intense 
absorption of LH2 at 850 nm (the so called B850 ring). 
Due to the relatively small distances between the chro- 
mophores in the B850 ring the interaction between Bchl 
molecules plays an important role in determining their 
spectroscopic and functional properties. A second ring of 
9 (or 8) weakly interacting Bchls is bound by the amino- 
terminal domain of LH2 (chromophore-chromophore dis- 
tance of about 2.1 nm) and is largely responsible for the 
absorption at 800 nm (the B800 ring) . Most of the spec- 
troscopic results for B850 band can be well explained in 
terms of the exciton model invoking available structural 
data [23L[2^. However, the single molecular fluorescence 
data and studies of the temperature dependence of the 
absorption spectra reveal that the conventional exciton 
model has to be modified due to the exciton interaction 
with the protein surrounding. The modified (dichoto- 
mous) exciton model was proposed to explain these dis- 
crepancies [25l - [27j . In contrast, the absorption and ex- 
citation dynamics in the B800 ring are usually consid- 
ered in terms of localized excitations of individual chro- 
mophores, due to the weakness of the intermolecular cou- 
pling. Despite of that, signatures of weak excitonic cou- 
pling have also been identified in the optical spectra of 
photosynthetic antennae of LH2 [7] , single molecular ex- 
citation spectra at low temperatures [28], and 2D spectra 
of its analogous LH3 complex Similar effects can be 
expected in other molecular aggregates with weakly cou- 
pled chromophores, DNA stacks, polymers, etc. 

In this paper we formulate and develop a theoretical 
model, which allows us to describe exciton dynamics in 
the system of weakly coupled molecular aggregates in 
terms of the reduced density matrix (RDM) in the basis 
of the site representation of the chromophores. The ex- 
citon coupling is treated as a perturbation and equations 
of motion (EM) for the RDM are derived in the second 
order. This enables us to retain main excitonic effects 
such as transition dipole moment redistribution, transi- 
tion energy shift and a characteristic excited state ab- 
sorption (ESA) shift, while simultaneously working with 
localized states. 

The paper is organized as follows. In the next section, 
we introduce the model system and its Hamiltonian. In 
Section nm we derive EM for a general weakly coupled ex- 
citonic complex, and in Section IIVI we discuss long time 
limit of these equations, leading to Forster type of relax- 
ation. We discuss calculations of linear absorption and 
two-dimensional photo-echo spectra in Section |Vl Model 
calculations are presented in Section IVll 



II. MODEL HAMILTONIAN 

Let us consider an aggregate of K two-level chro- 
mophores. The electronic ground state of such an ag- 
gregate can be described by a product state 




Figure 1: Illustration of the excitonic effect in a homodimer. 
Electronic excited states |1) and |2) with transition energy 
fiojeg localized on the individual molecules (depicted by their 
transition dipole vectors di and d2) form delocalized eigen- 
states |-f) and |— ) of the total electronic Hamiltonian. The 
result of excitonic interaction in a homodimer is splitting of 
transition energies by twice the resonance interaction J, and 
consequently an offset of ESA with respect to ground state 
absorption. The delocalized states have new transition dipole 
moments d+ and d_ corresponding to a sum and a difference 
of di and 6,2, respectively. 



\g) = \gi . ..gx) = |.gi)|52) ■ ■ • \9k}, (1) 

where the state \gi) is the ground state of the i—th 
monomeric chromophore. Excited states of the aggre- 
gate can be constructed using monomeric excited states 
\ei). Thus, we will represent single- and double-excitation 
states as 

\n) = \gig2---en---gK), (2) 

\N) = \{n,m)) = \gi...e„...em---gK), n< m. (3) 

Excited states with more than two excitation will not 
be considered as they can often be neglected in the first 
(linear) and third order spectroscopic experiments that 
we have in mind here. Unless stated otherwise, we will 
use capital letters to denote a double-excitation index, 
e.g. A = {o,,b), and lower case letters to denote one- 
excitation states. The diad (a, b) will be used to denote 
two-excitation states when the knowledge of the under- 
lying one-excitation states is required. Also sums over 
double-excitation states will be used or ^^^^ ^-j with 

the meaning of J^a^i T,b=a+v 

Individual chromophores are described by their 
ground- and excited state electronic energies ef and e^. 



nuclear potential energy surfaces of the ground- and ex- 
cited states Vf{Q) and V^^{Q) and the nuclear kinetic 
energy Ti{P). If we assume for a while that individual 
chromophores in the aggregate do not interact we arrive 
at the following Hamiltonian 



Hq — Hb + Hei , 



(4) 



where 



K 



Hei = e,\g){g\ + (e„ + (F„(Q) - \n){n\ 



K(K-l)/2 

E {^N + {VN{Q)-Vg{Q)))\N){Nl (5) 

N=l 



and 



HB = {T{P) + Vgm\9){9\ 
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B=(1,4) 



A=(1,2) 



J =J 
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n=l 



K(K~l)/2 

+ E {VN{Q)-{VN{Q)-Vg{Q)))\N){N\. (6) 
Here, we introduced 

K K 

T{P)=Y,T,{P), V,{Q)^Y.^!{Q), (7) 



1=1 



i=l 



Vn{Q) = E ^/(^?) + "^niQ). 



Vf{Q) + V,^Q) + VnQ). 



(8) 



The mean value of the difference between the potential 
energy in the state denoted by index a and potential 
energy in the electronic ground state was denoted as 

(K.(Q) ~ VgiQ)) = TrQ{{Vc.{Q) - VgmWe,}, (9) 

where a = n, N . The density operator Weq represents 
the equilibrium state of the nuclear DOF in the electronic 
ground state. 

Introducing excitonic interaction, the total Hamilto- 
nian reads 



Figure 2: Coupling between two two-excitation states A = 
(1,2) and B — (1,4). The coupling J24 enables transition 
between the two states, while the one-excitation state 1 is 
shared. 



Hj = E Jnm\n){m\ + E Jnm\N){M\. 



(11) 
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For two distinct double-excitation states M — (m, n) and 
N = (fc, I) we assume the following ansatz for the reso- 
nance coupling 



JmN — J(m,n)ik,l) — ^mkJnl + S'nlJmk 
^mlJnk ~t- ^7ikJral- 



(12) 



The ansatz is schematically described in Fig. [5] which 
demonstrates how two two-exciton states transfer from 
one to another by a transition of one excitation, while 
the other is shared by both double-excitation states. 

It is important to note that due to the resonance cou- 
pling Hamiltonian, iJj, neither the single- nor double- 
excitation states are eigenstates of the total electronic 
Hamiltonian. The actual eigenstates of the electronic 
Hamiltonian (obtained by its diagonalization) can be ex- 
pressed as linear combinations of the single- and double- 
excitation states, and they are usually termed one- and 
two-exciton states, respectively. 



III. EQUATIONS OF MOTION 
A. Projection operator method 



H — Ho + Hj, 



(10) 



We start investigation of the dynamics of the molecular 
aggregates with the Liouville-von Neuman equation for 
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the total density matrix W{t), 

^W{t) = -tCW{t). (13) 

Here, we defined the Liouville superoperator (or Liouvil- 
lian) £A = i [H, A] . The total Liouvillian can be written 
in terms of the system and the resonance interaction Li- 
ouvillians 

C = Co+Cj. (14) 

The solution of Eq. ([13]) can be written in terms of the 
evolution superoperator U{t) — exp{~jrCt} as W{t) — 
iJ{t)W{0). In interaction picture with respect to Cq, we 
define 

and 

W^'\t)=ul{t)W{0), (16) 

where Uo{t) = exp{ — ■|£oO- Using the usual projection 
operator method (P^ = P, Q = 1 — P) we get 

^PW^(^)(t) = ~tPCj{t)PW^'\t) - iPCj{t)QW{to) 
t 

- j dTPCj{t)QCj{T)PW^'\T), (17) 
to 

where we have already truncated the perturbation expan- 
sion in Cj at the second order. Although Eq. p7)) is valid 
for an arbitrary projector P, for the quality of the second 
order approximation, the choice of projection operator P 
is crucial. The best guidance to this choice is provided by 
the physical conditions at which the equation is applied. 
In optical spectroscopy, we often deal with systems that 
are in the electronic ground state at the initial time, and 
the bath DOF are relaxed into the canonical equilibrium. 
Thus, W{tQ) = \g){g\Weq = PggWeq- According to the 
Condon principle, the bath part of this initial condition 
is not changed upon an ultrafast photo-excitation and 
the initial condition for relaxation of the nuclear DOF is 
still given by a canonical density matrix. As a result, the 
projection operator 

PA^TrQ{A}W,q, (18) 

has the convenient property QM^(to) — 0, and the so- 
called initial term iP C j{t)QW {ti^) is identically equal to 
zero. This is equivalent to the statement that the pro- 
jection by P does not lead to a loss of information about 
the system at to- 

Another physically important situation occurs when 
the system has already spent some time in the electroni- 
cally excited state, and the bath DOF at an excited elec- 
tronic state I a) has relaxed into local equilibrium repre- 
sented by the density matrix W^^. The projection oper- 
ator which eliminates the initial term is 

PA = Y.(a\TrQ{A]\a)W:^\a){a\. (19) 

a 



This discussion shows that, in general, it would be desir- 
able that the projection operator P were time-dependent. 
It is indeed possible to formulate corresponding time- 
dependent projection operator technique rigorously (see 
e.g. Ref. ^]). In this contribution, we will not discuss 
such a time evolution of the projection operator, and we 
will keep in mind that the present formulation is valid 
only for not too long relaxation times. We will how- 
ever discuss Eq. ([T7]) with projection operator, Eq. (ITUl) . 
later in Section IIVI to demonstrate the relation between 
the theory developed here and the FRET. 

The main feature of the present theory is that the reso- 
nance coupling Hamiltonian Hj is the term with respect 
to which we apply perturbation theory. If the perturba- 
tion theory is applied with respect to the system-bath in- 
teraction Hamiltonian with projection operator given by 
Eq. (jlSp , we arrive at the Redfield type relaxation equa- 
tions (see e.g. [30]), while using the projector, Eq. (IT^ . 
leads to so-called Modified Redfield equations [31. ,32,]. 

B. Projector vs. interaction picture 

Eq. pT|) is an equation of motion for the RDM. How- 
ever, the projection operator is applied to the interaction 
picture density matrix W''^\t) and not to W{t) as one 
would expect. We therefore need to express the evolu- 
tion of p{t) = TrQ{W{t)} in terms of the evolution of 
pit) — TrQ{W'^^\t)}. Because Hamiltonian operators 
Hei and Hb act on different Hilbert spaces, they com- 
mute and we can write 

p{t) = ul{t)TrQ{ul{t)W{t)UB{t)}U,i{t). (20) 

Matrix elements Pab{t) of the RDM therefore read 

Pabit) = e^^-^'TrQ{Ul{t)Wabmb{t)}. (21) 

Due to the properties of the trace operation, we find that 
for populations 

Paa{t) = Paa{t). (22) 

For the coherences Pab{t), a ^ it is not possible to write 
directly such a simple result. However, an approximate 
relation between the two RDMs can be established by 
studying the case where Hj = 0. In this case, Wab{t) = 
UB{t)WeqUg(t) and pab{t) — Pab{0), becausc all the time 
evolution was accounted for by the interaction picture. 
At the same time, however, we can show that in second 
cumulant approximation 

Pabit) - Pa6(0)e-^"-*-»°(*)-»^(*), a ^ 6, (23) 

where ga (t) is the well-known line shape function associ- 
ated with the transition \ga) ]ea)l2y|- Thus, because 
the purpose of the interaction picture is to suppress the 
time evolution due to Hamiltonian Hq we can write 

Pabit) = e*--*+(i-^-)[s<"(*)+f?(*)lp„,(t), (24) 
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so that the exponential prefactor compensates the "J- 
free" evolution of the RDM elements in the spirit of the 
interaction picture. 



C. Reduced density matrix equations 



We can now rewrite Eq. (|17[) in second order in terms 
of p{t) as 



ip{t)TrQ{We,Hj{t)}W,g + R{J^), (25) 



where 



z 

Ri J^) = j dTPCj[t)QCj{T)PW^'\T). (26) 

to 

In this subsection, let us use the lower case indices to de- 
note both the one- and two-excitation states. The traces 
in Eq. ([^5]) can be easily evaluated in second order cu- 
niulant approximation 

{a\TrQ{H,j{t)W,q}\b) = Jaf,e^"-*rrQ{[/t(<)t/fc(i)W^e J 



= J„be''""''*-(^-^"')[»"(*)+^"'(*)l = Jab{t), (27) 

and we arrive at 
d 



dt 



Pab{t) = X! Jac{t)pcb{t) 



^^p,,(t)J,f,(t)-i?(j2)/W^,,. 



(28) 



In the second order term, Eq. (pS)) . we have to evaluate 
two commutators of Hamiltonian Hj with the reduced 
density matrix 



+TrQ{Hj{T)W,,}p{T)TrQ{W,,Hj{t)} 
+p{T)TrQ{W,qH.,(T)H.j{t)} 

- p(T)TrQ{W,,H.,[T)}TrQ{W,,Hj{t)}\ . (29) 
The matrix elements involved in Eq. read in detail 
{a\TrQ{Hj{t)Hj{T)We,}\h) = ^ Jac Jcte*"°^*+*"^'" 

C 

X TrQ{Ul{t)U,{t)Ul{T)Ub{T)We,}. (30) 
{a\TrQ{Hj{t)\c) . . . We^{d\Hj{T)}\h) ^ Jac ■ ■ ■ Jdb 



X e-"=*+---TrQ{C/t(i)C7,(t)VF,,t/](r)[/,(r)} (31) 



and 



(a\TrQ{WeqHj{T)Hjit)}\b) = ^ JacJcbe''^-^^ 



X TrQ{UliT)UrXT)UUt)Ubit)W,,}. (32) 
We introduce an auxiliary function 

Mabcd{t,T)=TrQ{Ul{t)Ub{t) 

X Ul{T)Ud{T)Weq}e'^''*+'^^^\ (33) 
with the property 

Mabcd{t,T) ^ MS,,,{T,t) (34) 

and after changing the integration variable to r' = t — r 
we can write the EM in the form 

d i 

-QlPo.b{i) ^ -^Yl Jac{t)Pcb{t) 



R{J^)/Weq ^^Jdr TrQ{H.j{t)H.j{T)W,q}p{T) 
to 



t-to 

^$^Pac(t)Jc6(t)-E^ I dr 



-TrQ{Hj{t)Weg}TrQ{Hj{T)Weg}p{T) 



[JacJcdMaccdit, t - t) - Jac{t)Jcd[t ~ t)] pdb{t - t) 



-TrQ{H.j{t)p{T)W,,Hj{T)} 



[JtaJXabdit. t-r)~ JLitmt - r)] Pcd{t - r) 



+TTQ{Hj{t)Weq}-p{T)TrQ{W,qHj{T)} 



[JdbJacMdbac(t, t~T) ~ Jdb{t)Jac{t - t)] Pcd{t - t) 



-TrQ{Hj{T)p{T)WeqHj{t)} 



+ [JbdJdcMbddc{t^ t-^) 
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JMit)J:cit-r)]pac{t-T). 



(35) 



This is a form suitable for introducing both long time 
limit (^0 — > oo) and the Markov approximation {pab{t — 
t) ~ pab{t))- We apply only the later one, because we are 
interested mostly in short times. Now, it only remains to 
evaluate the four-index matrix Mabcdit,T) which is done 
in Appendix A using the cumulant expansion. In Markov 
approximation, we can write our equations as 



d_ 
di 



Pab{t) 



Jac{t)Pcb{t) 



J2 Pac{t)Jcb{t) - J2 [Raccdit)pdbit) " Kabdi*) Pcd{t) 
cd 



Rdbac(t)pcd{t) + Rbddc{t)Pac(t) 



(36) 



where we introduced relaxation tensor 
t 



Rabcd{t) = 



dr 



JabJcdMabcd{t,t — t) 



Jab{t)Jcd{t - t) 



(37) 



We derived EM of the RDM in the interaction picture. 
We use Eq. ([M)) to transform the RDM in Schrodinger 
picture if necessary. The first two terms in Eq. ([M)) 
correspond to the delocalization effect of resonance cou- 
pling given by Eq. (P7)) . While the term e*"^"'* in Eq. 
([27| originates from the interaction picture with respect 
to the electronic Hamiltonian, the presence of the line- 
shape functions in Eq. (|27)) shows that the magnitude of 
this coupling decreases exponentially with growing time 
t. Thus the bath fluctuations dynamically destroys the 
resonance coupling. 



D. Two-excitation states 

In higher order spectroscopies, ESA contributes con- 
siderably to the signal. A delicate balance of ground 
state bleaching and stimulated emission on one hand, and 
the ESA on the other hand, is behind the disappearance 
of the 2D crosspeaks when resonance coupling goes to 
zero. For weakly coupled aggregates, correct description 
of ESA, and correspondingly the dephasing of coherences 
between one-excitation and two-excitation state is indis- 
pensable. 

In this subsection, let us again denote one-excitation 
and two- excitation states by the lower case and the upper 
case letters, respectively. The evolution of the system in a 
two-exciton state A = (a, h) is described by the evolution 
operator 



C/a (t) -cxp<j --H^i 



exp|--^[i7a(8)l{a} + l{fc}«'m]t|, (38) 



I.e. 



UA{t) =Ua{t)®Ub{t)®l{aM- 



(39) 



We denoted the direct product of unity operators from 
Hilbert spaces except of those in a set {a, &,...} by 
l{a,b,...}- Expressions containing a product U Ait)U'^B^t) , 
where one of the excited states is shared by the two- 
exciton states A and B = (a, c), thus yield 

U(a,b){t)ul^,^(t) = ® Ubmlit)- (40) 

Considering the first order term and a general case M — 
(to, n), N = (fc, I) one arrives at 

JMNit) = JMNe'"'^'"'TrQ{Ulit)UNit)Weg} 



= 5mkJnl{t) + SnlJmk{t) + 5mlJnk{t) + SnkJ-mlit)- (41) 

For the terms in the second order of J we have in a 
complete analogy 

JMNJcdMMNcd{t,T) = 



SmkJ7ilJcdMnlcd{t,T) + SnlJmkJcdMmkcd{t,T) 



+ 5mlJ7ikJcdMnkcd{t,T) + 5nkJmlJ cdMmlcd{t, t) , (42) 

and consequently 

RuNcdif) = SmkRnlcdif) + 5nlRmkcd(t) 



+ 5mlRnkcd{t) + SnkRmlcd{t)- 



(43) 



Thus, all quantities corresponding to the two-excitation 
states can be expressed directly using the one-excitation 
quantities. 



E. Homogeneous limit 

To evaluate the relaxation tensor we need to evaluate 
the following two expressions 

t 

R'abcdit) = -^JabJcd j dTMabcd{t,t-T), (44) 



and 



Kbcdit) = ~Y^Jab{t) I drJUt - r). (45) 
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In order to simplify the equations, we will assume so- 
called homogeneous limit, where we have 



ga{t) = TJ. 



(46) 



This simple formula allows us to evaluate all terms in the 
relaxation matrix analytically. First we observe that 



and therefore 



(47) 



1 



(Tc + r<j) - iuJcd 



If the dephasing constants T are non-zero, no such prob- 
lem can occur with Eq. (PSI) . 

We stress that the homogeneous limit is used here 
for demonstration purposes only. The line broadening 
function, Eq. pS)) . corresponds to a limit of ultra- 
fast stochastic bath with correlation function C{t) — 
Td{t). Such correlation function does not allow for in- 
troducing required thermodynamic properties C{uj) = 
g?iw/fcBT(^^_^-j ^YiQ corresponding spectral density 
C{uj) = J^^dtC{t)e^'^* . Consequently, the relaxation 
tensor given by Eqs. (^5)1 and ([M)) does not lead to any 
finite temperature thermal equilibrium. We will treat cal- 
culations with realistic correlation functions elsewhere. 



(48) 



The R' elements are obtained in a similar manner. We 
start with a splitting of the AI functions 



X M:^,,(r)e--*+^"-^ 



(49) 



Such a splitting is possible for an arbitrary g(t) function 
and is not limited to the homogeneous limit. In general, 
the integrals, Eqs. (|44|) and (j45|) . can be evaluated using 
the Fourier transform. In homogeneous limit we find that 



(50) 



where 



aabcd = (1 - Sac + Sad)Ta + (1 - She - Sbd)'rb, (51) 



Pabcd = {Sad - Sac)Ta + {Sbc - Sbd)Tb, (52) 



and 



labcd = Tc -f Frf -I- {Sac - Sad)'ra + {Sbd - Sbc)Tb- (53) 

Using the definition, Eq. (|44|) . we get 

^abcd\'^) — 'Jab'Jcd'^ 



1 



labcd — Pabcd — i^cd 



(54) 



The case when '^abcd ~ Pabcd ~ i'-^cd = Oi which can occur 
for homo aggregates, has to be considered separately. Ac- 
cording to Eqs. and (I49p . when the denominator is 
equal to zero, the dependence on the integration variable 
disappears and the integral leads to t. Consequently, 



IV. LONG-TIME LIMIT OF EQUATIONS OF 
MOTION 



In order to establish the relation between our EM, 
Eq. ([55]) . and standard EM used to describe dissipa- 
tive dynamics and energy transfer, we show that the 
above derivation leads to the well-know Forster resonance 
transfer rates in the long time limit. Considering the pro- 
jection operator, Eq. we first find that 



PCj{t)PW'-^\t) = 0, 
which leads also to 

PCj{t)QCj{T)PW^'\T) 

= pcj{t)Cj{T)pw^'HT). 



(56) 



(57) 



Using the definition of the projection operator, Eq. (1191) . 
Eq. (Ull) turns into 



^Paa(t) 



X / dr 

to 



{C*ba{t'T)+Cba{t-T)}paa{r) 



{Cab{t-T)+C:,{t-T)}pbb{T) , 



(58) 



where 



Cab{t) = Trg [Ua{t)Ul{t)W,\'] . (59) 

The initial state of the bath W^"i' is the one in which 
chromophore b is vibrationally relaxed in the electroni- 
cally excited state, and a is relaxed the electronic ground- 
state. Thus, W^^ = WlqW^q. We apply the second order 
cumulant expansion and get 



Kbcd{t) = J,6J,,te»(--+--)*-(— +'^--)*. (55) Gab{t) = TvQ {Ua{t)Ul{t)w:i^] Ttq {Ug{t)ul{t)w!^^} 
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g-9o(t)-«"agtg-sJ (t)+i(wbg-2At,)t 



(60) 



where we used the property 2Aa — TrgjAVaW^*} — 
TrQ{l\VaWl'i} (see e.g. Ref. [20]). By substitution 
t' = t — T ^ hmit tQ ^ —oo and Markov approximation 
Paa{t — t') K, paa{t) we obtain relaxation rate equation 
with Forster rates (see e.g. Ref. [s^ ) 



Ka^b = 2i?e 



+ ^ ^ Raddcit) Regit). 



(64) 



c d^a,c 



Also, when both the single and double excited states are 
present in the Rabcd matrix, many elements of the matrix 
are zero, most importantly those where e.g. a and b are 
a single excitation and double excitation indices, respec- 
tively. For the coherences involving the one exciton and 
two exciton states we have EM 



d_ 
di 



PaB{t) = Jac(<)PcB(t) + 't^. P '^^ [t) J C B {t) 



(61) - ^ Raccd{t)PdB(t) + [KaBoit) + RoBacit)] Pcvit) 

cD 



cd 



This demonstrates the relation of Eq. to the Forster 
energy transfer rates. It is important to point out that 
the only difference between the two sets of equations is 
in the "initial" condition set on the bath part of the den- 
sity matrix by the choice of the projection superoper- 
ator. Thanks to our choice of the projection operator, 
equations used in this work retain the description of the 
electronic coherences, while Forster rate equations give 
no prescription for them. 



- X! ^*BDDci^)Pacit)- 
CD 



(65) 



It is possible to rewrite this equation entirely using the 
one-excitation indices by considering Eqs. (HT|) and 
Setting B = {a,TT), D = (7, 6) and C = {a, (3) yields 



V. ABSORPTION AND TWO-DIMENSIONAL 
CORRELATION PHOTON ECHO SPECTRA 

To calculate optical spectra we concentrate on coher- 
ence elements of the RDM that correspond to optical 
transitions. We apply perturbation theory with respect 
to electric field of an incident light to calculate absorp- 
tion and two-dimensional correlation photon echo spectra 

M. 



A. Optical coherences 

To calculate response functions needed for evaluations 
of optical spectra in general, we need first to calculate 
evolution operators U{t) which fulfill the relation 



pit)=U{t)piO). 



(62) 



We use a simple consequence of this equation, namely 

Uabcd{t) = Pab(i), (63) 

where Pab{t) is calculated using EM, Eq. with initial 
condition /Oab(O) = SacSbd- 

If any of the indices a, 6, c, d equals g, we have Rabcd — 
0. For the optical coherences involving the ground state 
we therefore obtain the following EM 



d_ 

di 



Pag{t) 



- Raccd{t)pdia^) (t) +Tl+T2-%, (66) 



cd 



where last three terms 71, T2 and Ti are somewhat 
lengthy. We present the first term here, 

K 7T-1 

7l = 1^ ^ Pa(crP)it)Ji3Ti{t) + Pa(aiT){t)Jaa{i) 



13=17+1 



K 



+ X! Pa(i^l3){'t)Jfi'y{t) +Y Pa(aa){t)Ja7,{t) 



(67) 



and the remaining two are presented in full in the Ap- 
pendix. From the point of view of simulation feasibility, 
Eq. represents the main advantage of treating a 

weakly coupled excitonic systems in the site basis, as op- 
posed to the treatment in the excitonic basis. Although 
Eq. (|66p is rather lengthy, one is concerned only with 
tensor quantities with the number of elements propor- 
tional to N^, where N is the number of chromophores, 
as opposed to ~ iV^ which would be required in excitonic 
basis. 



B. Absorption spectrum 

The absorption spectrum is given by expression 



a{uj) 



n{uj) 



Re / dte" 
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x{^dgal^agbg{t)dbgPggj, (68) 

where (...) represents an averaging over isotropic distri- 
bution of orientations of the molecular transitions with 
respect to the light polarization. The transition dipole 
moments dag have to be understood as projections of the 
transition dipole moments on the light polarization vec- 
tor e, i.e. dag = dag ■ e. The averaging is done over a 
product of two of such quantities. We have 

^ab = {{dag ■ e)(dfcg • e)) orient. = ^ ' (^9) 

O IdagllClhgl 

If one now defines dag to dag = \<iag\ one can write 

oo 

a(w) « -^i?e / dte*"^* 
n(w) J 



X ^abdgadbgUagbg{i)Pgg- (70) 

ah 

We use Eq. [70] in subsequent simulations of absorp- 
tion spectra. It is important to note that because we 
do not work with electronic eigenstates one cannot as- 
sume the so-called secular approximation {Uabcd{t) ~ 
5acSbdUabab{t)) to be valid, and the orientational factor 
does not reduce to simple 1/3. 
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Figure 3: Illustration of excitonic splitting in a model homod- 
imer system. Linear absorption spectrum of a homodimer 
with perpendicular transition dipole moments calculated us- 
ing Eq. (|70ll . The following parameters were used to illus- 
trate the influence of resonance coupling: ei = £2 = 12500 
cm~^, resonance coupling J — 50, 100, 150 and 200 cm~^, 
and r = 1/400 fs"\ 



and 



C. Two-dimensional spectrum 

Two-dimensional (2D) Fourier transformed photon 
echo (FTPE) spectroscopy is well described by the third 
order time dependent perturbation theory with respect 
to light-matter interaction [23| . The spectroscopic signals 
are expressed in terms of response functions correspond- 
ing to light-matter interaction events. These response 
function allow us to calculate an arbitrary third order 
response of a multi-level electronic system, provided we 
know the evolution superoperators lA and the transition 
dipole moment elements dy , for all involved electronic 
levels. The response functions are directly proportional 
to the observed signal if the incident pulses are infinitely 
short. Expressions for all response functions Rig and Rif 
(i = 1, . . . 4) involved in the calculations of the 2D FTPE 
of our model systems are summarized in Appendix IVIII 

In this paper we will consider only 2D spectra with 
zero population time, t2 — 0. The impulsive limit signals 
in the rephasing and non-rephasing configuration can be 
obtained as follows (see e.g. [SSj) 

Suits, ti) — i?2g(i37 0, ii) 



+ i?4g(i3,0,ti) -i?2/(t3,0,ii). (72) 
The total signal is given by 

sit3,h) = eih)SRit3,h) + e{-h)SNRit3,~h), (73) 

and consequently, the 2D spectrum, which is defined as 
a double Fourier transform of the signal S{t3, ti) is given 
by Q 



OO CXD 







CO oo 



dh J dtiSNR{t3,ti)e 



ibJsta+ibJiti 



(74) 







Rsgih, 0, ti) — Rlf{t3,0, ti 



(71) 



In each response function component of the 2D spectrum. 
Fourier transform can be performed in ti and is times 
separately. 



VI. RESULTS AND DISCUSSION 

The theory developed in the above sections has been 
implemented in the spectroscopic package NOSE [U 
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Figure 4: Illustration of oscillator strength redistribution due 
to resonance coupling in a model heterodimer. The mutual 
position and orientation of the transition dipoles are described 
by angles cjii and (/>2 from Fig. [T] The parameters of the 
model are ei = 12500 cm~\ €2 = 12600 cm^\ J = 50 cm"\ 
r = 1/100 fs-\ (jii = 7r/2 and 02 = 0, 7r/4, 7r/2, 37r/2 and vr. 

which was used to perforin simulations of impulsive limit 
2D spectra of several small model systems. We have cho- 
sen systems where effects of weak excitonic coupling, such 
as those reported in Ref. [l3| for LH3 could be expected. 
We demonstrate below that our local basis description is 
sufficient to account for such effects. To take advantage 
of analytic equations derived in earlier sections, we stay 
in homogeneous limit. Simulations taking advantage of 
the full description, including a finite bath correlation 
time will be presented elsewhere. 

A. Molecular dimer 

The simplest system where a weak excitonic coupling 
effect in excited state absorption can be observed is a 
molecular dimer. Resonance interaction leads to the 
splitting of the excited states, redistribution of the tran- 
sition dipole moments and a shift of excited state absorp- 
tion. These effects will be demonstrated here. In addi- 
tion, one can also expect energy transfer between the two 
split excitonic levels, formation of a coherence between 
excitonic levels upon excitation by light and its dephas- 
ing. This class of effects is associated with the evolution 
of the system in the excited state band, and will be stud- 
ied within our model elsewhere. Fig. [T] illustrates the 
dimer geometry and its excitonic splitting. In all dimers 
considered here, transition dipole moments lie in a z = 
plane. Orientation of the dipoles with respect to x and y 
axes is determined by an angle such that dx = \d\ cos 0. 

Absorption spectrum of the dimer displays the split- 
ting of the levels, as well as the transition dipole mo- 
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Figure 5: Influence of dephasing on excitonic splitting in a 
model heterodimer. The parameters of the model are ei = 
12600 cm"\ £2 = 12500 cm~\ J = 50 cm"\ = 7r/2 and 
02 = TT. The dephasing rates are V = 1/100, 1/200, 1/300 
and 1/500 fs'V 



ment redistribution. Fig. [3] presents absorption spectra 
of a model homo-dimer with resonance coupling vary- 
ing from 50 cm~^ to 200 cm"-'^ and dephasing parame- 
ters r = 1/400 fs-i. Because the monomeric transition 
energies of the two levels are the same, excitonic mix- 
ing of the two levels is maximal at any resonance cou- 
pling value. We have artificially chosen the dipole mo- 
ments perpendicular to each other to eliminate the effect 
of transition dipole moment redistribution. We can see 
from Fig. Othat the prediction of the absorption maxima 
agrees rather well with the prediction of excitonic model 
(splitting of 2 J). It can also be noticed that the split- 
ting is smaller than predicted by excitonic theory when 
resonance coupling is small, most likely due to the bath 
suppressing the exciton coupling term in Eq. (j27p . 

The effect of the transition dipole moment redistribu- 
tion is illustrated on Fig. U) A hetero-dimer with dif- 
ference of 100 cm~^ between the transition energies on 
the two monomers was chosen, and the absorption spec- 
trum was calculated for a fixed resonance coupling value 
of 50 cm~^. Different mutual orientations of the dipole 
moments lead to enhancement of the absorption on one 
or the other split level, depending on mutual orientation 
of the molecules. 

In Fig. Owe demonstrate that increasing the dephasing 
rate F leads to broadening of the absorption spectrum. 
Unlike in case of exciton splitting were the position of 
the line does not shift, here we can observe a small shift 
towards less pronounced splitting with increasing the de- 
phasing rate. 

The effect of excited state absorption offset cannot be 
demonstrated on an ordinary absorption spectrum. We 
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Figure 6: 2D spectrum of a homo dimer. Subfigures: (A) 
in line configuration with excitonic coupling J = —80 cm~^, 
disorder width A — 100 cm~^; (B) sandwich configuration 
with excitonic coupling J — 40 cm"'^ and A = 100 cm~^. 
Full line contours represent positive values from 10% to 100% 
of the maximum. Zero and negative contours are dashed, and 
spaced by 10% of the positive maximum. 
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have therefore calculated 2D electronic spectra at popu- 
lation time ^2 = for two dimer configurations. An in 
line configuration, Fig. [SJ'V, corresponds to two transition 
dipole moments oriented head-to-tail with the distance 
and dipole moment length chosen such that the dipole- 
dipole coupling leads to J = —80 cm~^. This results 
in an offset of the ESA towards higher frequencies. The 
sandwich configuration with two parallel dipole moments 
and the same center to center distance results in positive 
coupling J — AO cm~^ . The ESA appears on the lower 
frequencies (Fig. [6l3) in this case. Both calculations are 
performed with a diagonal Gaussian disorder with the 
FWHM of A = 100 cm-i. 



B. Small aggregates 

We have also investigated trimers, tetramers and pen- 
tamers. As our work is motivated by highly symmetric 
homo aggregates like LH2 we calculated 2D spectra of 
aggregates of N monomers with circular N— fold sym- 
metry. Dipole moments are all in plane with the ring 
formed by the monomers, and we assume an angle a be- 
tween the tangent touching the circle at the position of 
the monomer and its transition dipole moment. We com- 
pare two cases: a — (tangential orientation of the chro- 
mophores) and a — —j (radial orientation with dipoles 
pointing towards to center of the ring). It was shown 
in Ref. [s^ that these two configurations have a dis- 
tinct position of the excited state absorption. Fig. [3^ 
presents 2D spectra of an average trimer with a — 0, 
Fig. [7j3 presents the same trimer calculated averaging 
over 100 realizations with energetic disorder of A = 100 
cm~^. The former figure reveals real part of the simple 
complex Lorentzian lineshape which is a consequence of 
the homogeneous limit assumed here. The excited state 
absorption is found below the ground state contribution 
in this case. For the radial configuration, i.e. a = 



Figure 7: Two-dimensional correlation photon echo spec- 
trum of a model trimer. Subfigures: (A) radial configuration 
a = — ^ , single realization; (B) radial configuration, averaged 
over disorder with A = 100 cm~^; (C) tangential configura- 
tion a — 0, single realization; (D) tangential configuration, 
averaged over disorder with A — 100 cm~^. All site energies 
are e — 12500 cm" and r = 1/300 fs"\ Contours as in Fig. 

El 

show on Figs. [7p and[7|D the excited state absorption is 
found above the ground state contribution. If a trimer is 
considered a member of the family of TV— fold symmet- 
ric aggregates this result is the opposite of the expected 
effect identified in Ref. [sl]. However, trimer has to be 
considered a special case with respective angles between 
the chromophores very different from the larger aggre- 
gates of the same symmetry. In larger aggregate we can 
expect that the ESA will be in a position similar to the 
in line dimer for tangential orientation, and in a posi- 
tion similar to the sandwich dimer for radial orientation. 
Indeed, already the pentamer follows the rule found for 
larger circular aggregates. As we can see on Fig. [5] the 
two configurations have now position of the ESA in agree- 
ment with Ref. [35]. 

Figs. [6] to [8] demonstrate that the theory developed 
in this paper reproduces correctly the ESA features of 
2D spectra of symmetric weakly coupled excitonic aggre- 
gates. From the theoretical point of view ESA features 
are a result of a delicate balance between ESA and GSA 
contributions which cancel exactly in case of uncoupled 
chromophores. This feature makes 2D of uncoupled chro- 
mophores additive. 



C. Additivity of the 2D spectra 

We will demonstrate that our theory fulfills the addi- 
tivity property, and show that the GSA and ESA contri- 
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Figure 8: Two-dimensional correlation photon echo spectrum 
of a model pentamer. Subfigures: (A) radial configuration 
a = — ^ , single realization; (B) radial configuration, averaged 
over disorder with A — 100 cm~^; (C) tangential configura- 
tion a = 0, single realization; (D) tangential configuration, 
averaged over disorder with A — 100 cm~^. All site energies 
are e — 12500 cm" and r = 1/300 fs"\ Contours as in Fig. 

El 



Figure 9: Demonstration of the additivity of 2D spectra. Sub- 
figures: (A) Full 2D spectrum of a tetramer consisting of four 
uncoupled molecules, (B) the same calculation with ESA ig- 
nored, (C) Tetramer with only the lowest and highest energy 
molecules coupled to each other by J = 100 cm~^, (D) The 
same case with ESA ignored. Site energies are e — 12800, 
12600, 12400 and 12200 cm"\ disorder width A = cm"^ 
and r = 1/150 fs~^. Contours as in Fig. El 



buttons are correctly balanced. To this end we consider 
a hetero tetramer composed of four chromophores with 
distinct transition energies oji = 12600 cm"^,CL>2 = 12500 
cm-i,a;3 = 12400 cm^^ and uj^ = 12300 cm-\ Fig. 
presents a 2D spectrum of uncoupled tetramer. No cross- 
peaks and negative features appear and the spectrum 
is a sum of monomeric 2D spectra. When only ground 
state to one-exciton band transitions are considered and 
the aggregate ESA is ignored many crosspeaks appear in 
the 2D spectrum (Fig. EP)- AH these crosspeaks are ex- 
actly canceled by the ESA contribution. When the low- 
est and highest energy monomers are coupled (here with 
J — 200cm^^) the 2D spectrum becomes a sum of two 
independent monomers and a coupled dimer. Fig. |9p. 
Again, when the ESA contribution is removed, the redis- 
tribution of cross-peak amplitudes in ground state contri- 
bution can be clearly seen, but all cross-peaks involving 
the coupled dimer and the independent monomers are 
canceled out in Fig. [Sp. 

To calculate 2D spectrum of uncoupled monomers one 
can therefore simply calculate individual 2D spectra and 
sum them. However, the balance of ESA and the ground 
state contributions can be disrupted also by relaxation 
processes like FRET. Two monomers or two composed 
systems that are coupled weakly so that cross-peaks due 
to mutual interactions are too weak to be resolved, but 
nevertheless strong enough to enable energy transfer will 
show energy relaxation crosspeaks. This is e.g. the case 
of LH3 (see Ref. [14]). In case of LH3 one has to ac- 



count for additional stimulated emission from the states 
populated by energy relaxation, and other similar pro- 
cesses. All this is included in the EM presented in this 
paper. Moreover, not all states involved in the theory 
have to be localized on individual monomers. We can 
equally divide the system into parts where excitonic in- 
teraction dominates, and start from excitons formed by 
this interaction. Mutual interaction of such blocks is then 
described by the theory presented here. 



VII. CONCLUSIONS 

In this paper we have derived EM for the reduced den- 
sity matrix of a system of weakly coupled chromophores 
interacting with an environment. The weak excitonic 
coupling is treated in the second order perturbation 
theory and the environmental degrees of freedom are 
described within the second cumulant approximation, 
which for some type of systems provides an exact so- 
lution. We show that our equations are related to the 
Forster relaxation rates. In contrast to the usual Forster 
type equations, we provide a detailed prescription for the 
evolution of coherences. Thus, we are able to describe an 
effect of dynamic localization, where the bath destroys 
not only a wavepacket created in the complex by ultrafast 
excitation, but also the coherence established by the weak 
resonance coupling. By simulations of model systems in 
homogeneous limit we demonstrate that 2D spectroscopy 
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reveals the excitonic coupling by an offset of the excited 
state absorption, and that our local bases description of 
the systems dynamics is fully sufficient to account for this 
effect. 
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1 - iUbit) - Qi''\t)\ \i + iHc{t) ~ g[r\T) 



Tr, 



){l + naimbit) - naimdr) + UaimdiT) 



b{t)Uc{T) - Ubimdir) + nc{T)Ud{r)] 



Appendix A: Cumulant expansion evaluation of the 
Mabcd{t,T) matrix 

Expanding each evolution operator Ua{t) up to the sec- 
ond order, 



where 



Ua{t)^i-ina{t)-g^^\t), 
ul{t)^i + i'Ha{t)-gV{t), 



(75) 



(76) 



(77) 



- TrQ[g'^-\t) + g^^\t) + (?(-)(r) + g^;\T)]. (84) 

Assuming that gab{t) — if a ^ b , and taking into 
account that a ^ b and c ^ d in the M function we have 



Mabcd{t,T) = (*■^)+''^-*+*'^-^ 



where 



Fabcd{t,T) = ~gl{t) - gbit) ~ gl{T) ~ gair) 



-Sac {9a{t) - ga{t - t) + g*(r)) 



-Sad {9a{t) - ga[t ~ t) + gl{T)) 



(85) 



and 



z r 

^^^^^^^ ^hj'^^j dr'^V{T)/^V{T'), (78) 







Z T 



(79) 



the auxiliary matrix M can be evaluated in the second 
order cummulant approximation. To that end we evalu- 
ate the following traces 



TrQ{Ha{t)W,q} = Q, 



TrQ{gi+\t)W,,} = gait). 



(80) 



(81) 



+Sbc {9b{t) ~ 9b{t - t) + glir)) 



Sbd{gb[t)'gb{t-T)+gl(T)). (86) 



Appendix B: Third order response functions 

Here we present third order response functions for a 
three band system in state representation, using Einstein 
summation convention. The upper indices denote the 
bands so that g corresponds to the ground state band (for 
excitonic system only one state is assumed to be in the 
ground state band), e and / represent the one-exciton 
and two-exciton bands, respectively. Lower indices de- 
note the states within the bands. In all equations below, 
index a represents a ground state index and consequently 
a = g and pa = Pg, where pg is the initial population of 
the ground state. 



TrQ{gi-\t)W,,} ^ g:{t), 



(82) 
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X Ul^%{h)U%t{hn2'a{tl)Pa. (87) 



(eg) / 



K 

' E 

S=a+1 



if 

E 

5=7r+l 



ir-l 



E [R*caa^ii) + ^7-c(i)] Pc(7^) (<) • (95) 



7=1 



X^J;)(t3)<j,(t2)Wr,(^l)Pa, 
xZ^(;;/)(t3)<i,(i2)WL^:l(<l)Pa, (90) 

i?3.(i3,t2,ii)=(Krvi;''Virvif)) 



(89) 



K K 

E E 

a=l I3=a+l 



+ R*7racrl3 (*) + KpTra (0 ) Pa(a/3) (*) 



K K 



E E K55^(t)Pa(.^)(i) 



;3=cr+l S=a+1 



xZ^5,^](<3)<i(t2)Z^,ti(tl)p., 



i?3/(i3,i2,ti)=(i^,rv,(7Virvr)) 



xwg2(i3)z^}S(*2)wL^:l(ii)p., 



(91) 



(92) 



cr-1 if 

E E 

7r66a 

a=l S=<7+1 



K TT-l 

E E 

(i)Pa(7r/3) (0 

3=7r+l 7=1 



i?4,(i3,t2,ii) = (v^,rVi7Virv,(;^') 

xZ^S(^3)Z^}S(i2)^if (tl)Pa, (93) 

i?4/(i3,t2,ti) = (^,rVifVi/^)yi;^') 

xZ^(;^)(t3)Z^jS(*2)W(^gI(il)Pa. (94) 

The sign (...) represents orientational averaging over 
possible orientations of a molecular system with respect 
to the polarization axis of the incident light. The orienta- 
tional averaging is preformed for an isotropic distribution 
of orientations according to Refs. (sgI. [37j. 

Two-excitation Terms in Equation of Motion 



E E ^*'yita{i)Pa(a7T){i) 



a— 1 7—1 



+ E E -^^77" {t)Pa{aa) (t) 



a— 1 7—1 



K cr-1 

■ E E 

{t)Pa{^P)it) 

j3=(j+l 7=1 



K K 
■ E E -^^<5<5/3(OPa(7r/3)(i) 



In this appendix we present the details of 72 and Ts 
terms of the EM, Eq. (j66|) . for two-excitation states. 



K (T-l 



c=l 7=1 



E E ^^A'Aa (i)Pa(a7r) (0 • 



S=-K+l a=l 



(96) 
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